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Abstract 

Background: In radiation protection, biokinetic models for zirconium processing are of crucial importance in dose 
estimation and further risk analysis for humans exposed to this radioactive substance. They provide limiting values of 
detrimental effects and build the basis for applications in internal dosimetry, the prediction for radioactive zirconium 
retention in various organs as well as retrospective dosimetry. Multi-compartmental models are the tool of choice for 
simulating the processing of zirconium. Although easily interpretable, determining the exact compartment structure 
and interaction mechanisms is generally daunting. In the context of observing the dynamics of multiple 
compartments, Bayesian methods provide efficient tools for model inference and selection. 

Results: We are the first to apply a Markov chain Monte Carlo approach to compute Bayes factors for the evaluation 
of two competing models for zirconium processing in the human body after ingestion. Based on in vivo 
measurements of human plasma and urine levels we were able to show that a recently published model is superior to 
the standard model of the International Commission on Radiological Protection. The Bayes factors were estimated by 
means of the numerically stable thermodynamic integration in combination with a recently developed copula-based 
Metropolis-Hastings sampler. 

Conclusions: In contrast to the standard model the novel model predicts lower accretion of zirconium in bones. This 
results in lower levels of noxious doses for exposed individuals. Moreover, the Bayesian approach allows for 
retrospective dose assessment, including credible intervals for the initially ingested zirconium, in a significantly more 
reliable fashion than previously possible. All methods presented here are readily applicable to many modeling tasks in 
systems biology. 

Keywords: Bayesian inference, Model selection, MCMC sampling, Compartmental model, Internal dosimetry, 
Systems biology 



Background 

Radioactive zirconium (Zr) isotopes are produced in large 
quantities in nuclear fission reactors; one of the most com- 
mon fragments in uranium fission is the radioactive 95 Zr. 
For example, the estimated released 95 Zr activity of the 
Fukushima and Chernobyl accidents is considered to have 
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detrimental health effects not only via irradiation, but also 
via the intake of edibles [1,2]. The estimation of radiation 
doses is indispensable for risk analysis. This is true for 
occupational exposure [3] and patients undergoing diag- 
nostic and therapeutic nuclear medicine [4] as well as for 
the public in general [5]. To calculate the radiation dose a 
mathematical model is required for quantifying the depo- 
sition of radioactivity from the incorporated radionu- 
clide inside the human body. In internal dosimetry, this 
model is called biokinetic model as denned by the Inter- 
national Commission on Radiological Protection (ICRP) 
in [5]. Also, the ICRP put forward the current standard 
model, which we will simply denote the ICRP model. The 
parameters of this model were mostly derived from ani- 
mal data. In order to obtain more reliable dose estimates 
for humans, the Helmholtz Zentrum Munchen (HMGU) 
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developed a new, physiologically more plausible biokinetic 
model [6]. It is based on the processing of non-radioactive 
Zr isotopes in 16 investigations with 12 healthy human 
subjects. The measurements were taken in vivo in plasma 
and urine up to 100 days after administration by appli- 
cation of the double tracer technique. Moreover, a global 
statistical analysis method has been developed and the 
uncertainty and sensitivity of the HMGU model parame- 
ters were analyzed [7,8]. 

The biokinetic models mentioned above incorpo- 
rate basic processes in the human physiological system 
[3,5,9,10]. Mathematically, this is characterized as follows: 
All major human organs and tissues are simplified in the 
model structure as separate compartments that represent 
kinetically homogeneous amounts of radionuclides; the 
connections between organs and tissues are described via 
transfer rates, i.e. model parameters that represent the 
exchange rates between the individual mutually exclu- 
sive compartments. These multi-compartmental systems 
along with their transfer parameters describing the kinetic 
behavior of radionuclides in the human body are called 
compartmental models [5,11]. Throughout this paper we 
use the terms biokinetic model and compartmental model 
interchangeably. The transfer of substances into and out 
of compartments is governed by the law of mass balance. 
Transfer parameters are frequently evaluated on the basis 
of experimental data obtained from laboratory animals 
and, to a lesser extent, human beings [10]. Although ani- 
mal data is not directly comparable to human data, the 
former can nevertheless be very helpful as prior informa- 
tion. 

In this publication, we address the problem of model 
inference and model selection. A Bayesian approach 
enables us to cover model and measurement uncertainties 
for a compartmental model based on human data, while 
simultaneously taking into account the prior information. 
The Bayesian framework provides a fully probabilistic 
approach [12]. It is grounded on the probability distribu- 
tion of a problem specific parameter space conditioned on 
the given data. This specifies a measure of belief for all 
possible parameter values. Similarly - albeit not identical 
- to confidence intervals, Bayesian analyses provide cred- 
ible sets for the parameters at stake, holding regions and 
limits of high parameter probability [13]. However, as they 
are intrinsically driven by prior informations, some care 
has to be taken in their interpretation [14]. 

For an overall assessment of the two competing bioki- 
netic models for Zr, the previous model parameter uncer- 
tainty analysis [7,8] is not sufficient, because uncertainties 
arising from the model structure were not taken into 
account. This shortcoming was addressed by our Bayesian 
approach. Considering the models themselves as a ran- 
dom variable allows to compute the probability for each 
of the models conditioned on given data. The ratio of the 



marginal likelihoods of two models, i.e. the ratio of the 
probability for the data to be produced by the specific 
model, is known as the Bayes factor, a quantity intro- 
duced by Jeffreys [15]. Performing model selection using 
Bayes factors is particularly useful when dealing with a few 
models only. While classical model selection approaches 
using statistics such as the AIC or likelihood ratio tests are 
based on single best parameter estimates [16], the Bayes 
factor takes into account all possible parameters values 
and thus deals with model uncertainty and avoids over- 
fitting issues [17,18]. Moreover, in contrast to classical 
tests, the Bayes factor provides evidence for either one 
of the (possibly non-nested) models by definition. With 
the introduction of Markov chain Monte Carlo (MCMC) 
methods [19-21] as tools for sampling from probability 
distributions, a remarkable increase in Bayesian analyses 
was noticed. However, due to very complex probability 
surfaces these approaches often struggle with sampling 
efficiency [22]. In order to avoid resulting convergence 
issues of the MCMC approach, we combined a technique 
called thermodynamic integration with a novel copula- 
based Metropolis-Hastings sampler [23]. This provides 
numerically stable results for the inference process. 

The HMGU and ICRP models were compared based 
on in vivo plasma and urine data of 16 investigations of 
12 human subjects [6] using Bayes factors. More pre- 
cisely, the models were evaluated for each investigation 
individually and for the concatenated data of all investi- 
gations. The latter allows to infer transfer rates (including 
credible intervals) for an average individual. We further- 
more provide an analysis based on the (i) plasma data and 
(ii) urine data individually. Throughout the analysis, the 
HMGU model turned out to be superior compared to the 
ICRP model with respect to covering the specific data. 
This means the HMGU model more precisely explains 
the given observations and therefore the processing of zir- 
conium in the human body. We then used the HMGU 
model to analyze the accretion of zirconium in bones: 
not only did we observe a delayed aggregation but also to 
lesser extents than predicted by the ICRP model. Lastly, 
the Bayesian framework yielded credible bounds for ret- 
rospective dose assessment of an average individual, this 
is, based on the concatenated data of all 16 investigations. 
We provide a user-friendly estimation table for the predic- 
tion of initially ingested zirconium concentrations for ex 
post measurements. This impacts the estimation of intake 
and radiation dose, especially the bone dose, when aiming 
for personalized occupational monitoring data. 

Methods 

Biokinetic models for zirconium processing 

We infer the biokinetics of zirconium as it is processed 
in the human body. The currently used compartmen- 
tal model was recommended by the ICRP in [5,10,24] 
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(Figure 1A). It consists of eleven compartments and 15 
transfer rates. Zirconium enters the body via the stomach 
compartment 3/9 and is processed until it reaches any of 
the two final compartments urine, 3/7, or feces, y$. Some 
of the transfer rates and compartments of the ICRP model 
are however physiologically questionable: The direct mass 
transport from the two bone compartments to the urinary 
bladder contents and upper large intestine compartments 
or the distinction between trabecular bone surface and 
cortical bone surface as such. In order to address these 
shortcomings, we have recently proposed the alternative 
HMGU model [6] combining the two bone compart- 
ments into one single compartment and replacing these 
mass flows by physiologically more plausible transfer rates 
(Figure IB). Altogether both models share eight transfer 
rates, which we denoted by x\, ... ,x%. Transfers present 
in just one of the models have a unique index to facilitate 
distinction. 

The dynamics of both models are described by a sys- 
tem of coupled linear first-order ordinary differential 
equations (ODEs): For each compartment yj with time- 
dependent concentration yj(t), the rate of change of yj is 
given by 

aeAt fie Ay. 



where A^. is the set of indices of all transfer rates x a 
flowing into compartment yj, A~ the set of indices of all 
transfer rates flowing out of compartment yj and y[ Xa \ is 
the compartment which is connected to yj by the trans- 
fer rate x a . For instance A^ 5 = {8, 10}, y[ X8 ] = yio and 
y [xm] = y x in the HMGU model. We have y 9 (0) = 100% 
and therefore 3/^9(0) = 0%, this is, the complete amount 
of zirconium is initially contained in the stomach com- 
partment. The explicit model specific ODE systems can be 
found in the Additional file 1 sections 1.1 and 1.2. 

Experimental data 

The human biokinetic data was measured in a sta- 
ble tracer study executed at the Helmholtz Zentrum 
Munchen (HMGU) [6,25]. It includes 16 investigations 
of 12 healthy humans with ingestion of a investigation- 
specific amount of isotopically enriched stable zirconium. 
The administered amount was based on the individuals 
weight, aiming at a dose of 0.09mg stable tracer per kg 
body weight. Tracer concentrations were determined in 
blood plasma and urine. For the plasma data, samples 
were taken several times during the first day in increasing 
intervals, and more scarcely later on. Urine was collected 
completely in 12-24h periods on several days. The last 
samples were taken at lOOd after tracer administration. 
Tracer concentrations were normalized to the respective 
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Figure 1 Models for the biokinetics of zirconium. A: ICRP model. The model consists of eleven compartmentsyi,...,yn and 15 time independent 
transfer rates xi,...^8Ai3#---Ai9- B: HMGU model. The model consists often compartments yi ,...,yio and twelve transfer rates x-\ ,...,xn- In both 
models zirconium enters the body in the stomach compartment yg and diffuses through the system until it reaches either one of the two final 
compartments urine, yj, or feces, yg. The gray compartments yi and yi are directly related to the datasets measured. 
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tracer amount ingested in each investigation, such that 
the total ingested amount corresponds to 100% at t = 0 
in the stomach compartment 3/9. For model development, 
the transfer compartment was taken to be identical with 
blood plasma, and concentrations therein were expressed 
as % per kg plasma. The plasma concentration was scaled 
by the total amount of plasma in the body to get absolute 
concentrations [26]. Urine data was expressed as excretion 
rate in % per day. 

Model likelihood 

Technical limitations, such as missing in vivo measure- 
ment techniques for all involved compartments as well 
as noisy data introduce model uncertainties to biolog- 
ical systems [27]. Comparing different kinds of models 
based on single parameter estimates as done in maximum- 
likelihood approaches is thus very critical. For statistical 
model evaluation we here applied a Bayesian approach, 
taking into account the full parameter distribution: For 
each investigation i we assume that the data 



l > n i >i,l -i,2 - l > n i\ 



follows the solution c x k(t) of the differential equation 
given in (1) for any of the two models A4 k and some corre- 
sponding parameter vector x k , where the model index k e 
{H,I}. Here, M 1 is the ICRP model and M H the HMGU 
model. Corresponding to the notation in Figure 1A and 
IB, x 1 = (x\, . . .,x$,xi3, . . . ,#19) and of 1 = (%\, . . . ,#12). 
While y\ indicates measurements in plasma, i.e. in the 
transfer compartment y\ , y l j' designates measurements of 
the excretion rate in the urine compartment 3/7. There 
are n\ measurements in plasma and measurements 
in urine for investigation L Assuming furthermore that 
all data points contain normally distributed measurement 
errors, the investigation i and model k specific likelihood 
function has the form 



n i 

Ci(X k ,k\Vi) = fl 0^f|^fe),orA 



£?(***|Z>i) 



where c^(t a ) denotes the solution to Equation (1) for the 
transfer compartment y\ at time point t a , corresponding 
to the measurement at y l { a , for the parameter vector x k . 
Furthermore, ^c™ k (tp) is the derivative of the solution for 
the urine compartment 3/7 at time point tp, corresponding 



to the measurement yj , while o) is the univari- 

ate probability density function of the normal distribution 
with mean \i and standard deviation a. 

The standard deviations for plasma, of, and for 
urine, of, are fitted investigation-specifically by simulated 
annealing [28] before starting the MCMC sampling pro- 
cess. They correspond to the combined strength of all 
deviations from the ODE, e.g. to the size of the measure- 
ment error as well as to the magnitude of the difference 
between the ODE solution and the data points that is due 
to natural internal fluctuations not considered by an ODE 
approach. The biological variability between the individ- 
ual investigations is accounted for by the fact that we fit 
different of and of for each investigation i and thus get 
investigation-specific likelihoods. This leads to individual 
credible intervals for each parameter in each investigation 
in the MCMC sampling procedure later on. 

The complete data (i.e. concatenated data) likelihood is 
given by k\V) = Y\]ti k\V t ) for the com- 

plete data V = {£>i,..., £>i 6 } where in all Ci(x k ,k\Vi) 
the same fitted investigation independent of = o b and 
of — o u are used. 

For the calculation of £i(x k ,k\Vi) Equation (1) has to 
be solved based on x k . Since (1) is of first order, it can be 
written as 



dy x k(t) 

dt 



= A(x K ).y x ,(t), 



(2) 



where y x k(t) is the vector of all the compartments of 
model k and the time independent matrix A(x k ) repre- 
sents all the interactions between these compartments, 
depending on the transfer rate values x k . The analytical 
solution is then given by 

y xk (t)=e A ^ t . Yxk (t = 0). (3) 

The matrix exponential g^**)* was computed by eigen- 
value decomposition using MATLABs eig function (see 
Additional file 1 section 1.3). 

Bayes factors 

Specifying model specific, investigation indepen- 
dent prior distributions p(x k \k) based on combined 
human/animal data yields the posterior distributions of 
the parameters for model k according to Bayes' theorem: 

C i (x k ,k\V i )p(x k \k) 



p(x K \V h k) = 



(4) 



p(Pi\k) 

where p(T>i\k) denotes the marginal density obtained by 
integrating over the according parameter space X* Le. 

p(Vi\k)= f C i (x k ,k\V i )p(x k \k)dx k . (5) 
Jx k 

The quantity p(T>i\k) is called data evidence and is the 
basis for the model selection process. For comparing two 
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models k and k ', we can view the model index as a ran- 
dom variable, apply Bayes' theorem again and take ratios 
to cancel the denominator, yielding 



p(k\Vj) = p(Vj\k) p(k) 
p(k'\Vt) p{Vi\k>) ' p{k f ) 



(6) 



The ratio of the two likelihoods in Equation (6) is known 
as the Bayes factor 



B< 



k,k' 



P<Pi\k) 
P(Vi\k'Y 



(7) 



We had no initial preference for any of the models and 
thus chose a uniform model prior. The Bayes factor in 
Equation (7) then coincides with the posterior odds ratio 
between the models k and k f conditioned on the data T>i 
[18,29]. 

The Bayes factor compares the posterior probability 
p(k\T>i) that the data V{ have arisen according to model k 
versus the probability^// \Dj) = 1 — p(k\T>i) that V{ have 
arisen according to model k f , in contrast to single param- 
eter measures such as the AIC or the likelihood ratio test 
[16]. The models may be non-nested. Due to the eval- 
uation of each model on its whole parameter space X k 
(cf. Equation (5)), the Bayes factor naturally penalizes 
model complexity and thus prevents overfitting issues [30- 
32] . An interpretation of the likelihood ratio in the Bayes 
factor was given by Jeffreys [15], which can also be found 
in the Additional file 1 section 3. Most importantly, a 
Bayes factor B l k k , > 3 substantially favors model /<, while 
B l kk , > 100 decisively favors model /<. Clearly, for B l kk , < 1, 
model k' is favored with evidence B l k , k = 1/B l kk ,. 

Prior information 

Since the problem of radiation protection has been exten- 
sively researched, quite a large number of animal stud- 
ies have been conducted. These yielded excessive prior 
knowledge for the Bayesian modeling approach. As the 
ICRP model is the recommended model used for dose 
estimation, there naturally exists information on the dis- 
tribution types of the parameters involved in the model 
together with confidence intervals [7]. Four parameters 
have a lognormal distribution, five a triangular distribu- 
tion and six have a normal distribution (see Additional 
file 1 section 2.3 for details). From the confidence inter- 
vals, the parameters of the normal and lognormal distri- 
butions were calculated. For the HMGU model detailed 
prior information is also available from previous studies 
[7,8]. Here, eight parameters have a lognormal distribu- 
tion and four a triangular one (see Additional file 1 section 
2.3 for details). It is noteworthy that of the eight parame- 
ters shared in both models, x$ is the only one having dif- 
ferent distributions in the ICRP and HMGU model. Due 
to a lack of prior knowledge of the dependency structure 



between the parameters, the multivariate prior distribu- 
tion p(x k \k) of model k was taken to be the product of 
the univariate prior distributions p(x k \k) for each param- 
eter x k , i.e. p(x k \k) = \\ r p{x k \k). Each univariate prior 
distribution was truncated at zero. While Bayes factors 
were computed inter alia for each investigation separately 
(see Results and discussion), the same prior information 
was applied throughout all investigations. This is reason- 
able since the priors contain information from a huge 
variety of preceding experiments. 

Thermodynamic integration 

Computing the marginal of Equation (5) by numerical 
integration generally turns out to be a daunting task. 
There exist a variety of approaches to tackle this problem. 
Popular choices include Chibs method, which is based 
on a Gibbs sampling scheme and therefore not always 
easily applicable [33] or the Reversible Jump MCMC algo- 
rithm proposed by Green [34], based on an across model 
search. Since the latter involves sampling of an additional 
model parameter, the sampling generally takes longer than 
sampling within the different models only. We therefore 
applied a within model sampling technique called ther- 
modynamic integration (TI). It was proposed by Lartillot 
and Philippe [35] and Friel and Pettitt [36] and success- 
fully applied e.g. in Xu et al. [37]. The method splits apart 
the computation into several intermediate steps by intro- 
ducing an auxiliary "temperature" parameter r e[0, 1] 
that governs the influence of the parameter likelihood. 
The basis of this approach is the power posterior, which 
is the usual posterior modified such that the likelihood 
is exponentiated by the temperature parameter and then 
normalized accordingly to obtain a probability density: 



Px<Pi\k) 



Ci^j^ViYp^k) 

Z(Vi\k,T) 



(8) 



More precisely, the quantity of interest is the normaliza- 
tion constant 



z( Vi\k t r)= f C i (x k ) k\V i ) T p(x k \k)dx k 
Jx k 



(9) 



since it yields a way to compute the terms of the Bayes 
factor (cf. Equation (7)) by differentiating its logarithm 



d 
dx 



logz(Vi\k,r)= f log jCi(**,k\Vi) 
Jx k 



z{Vi\k,x) 
= E Pr [log £,(**, (10) 

and then integrating both sides with respect to r 

]ag(p(Pi\k)) = jf * E Pz [log£/(**»*|Z><)] dr, (11) 
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according to Calderhead and Girolami [38]. This means 
that the natural logarithm of the marginal likelihood 
can be computed as the integral over the expectation of 
the logarithmized data likelihood within the model with 
respect to the power posterior. The parameter r governs 
the flatness of the power posterior surface and, much as 
in the concept of path sampling [39], stabilizes the com- 
putation of Equation (5) [36]: choosing a discretization 
0 = ri < T2 < . . . < rjv-i < ?n = 1> we can compute the 
natural logarithm of the marginal likelihood p(T*i\k) by 
numerically approximating the integral in Equation (11) 

by 

N-l , 

\og(p{Vi\k)) « J2 5( T «+1 " T n)U Prn+l [log£/(**,*|Z>/)] 
n=l I 

+E Prn [log £i(x k ,k\Vi)] J. 

(12) 

Also, the expectation with respect to the power poste- 
rior in Equation (12) is approximated in the usual way by 
substituting it with the Monte Carlo estimate 

1 M 

E Prn [log dixKkWi)] log U (x\ my k\V^ 

m—l 

(13) 

where x^ denotes a sample drawn from p Tn (Di\k). For 
all our applications we chose a temperature schedule with 

N = 30 steps according to z n = (jS^ri) > n = h • >N 
to estimate \o%(p(T>i\k)) for each k and i as suggested by 
Calderhead and Girolami [38]. 

Copula-based Monte Carlo sampling 

The model, investigation, and temperature specific under- 
lying Markov chain Monte Carlo (MCMC) samples 
were drawn using the recently introduced copula-based 
Metropolis-Hastings (MH) algorithm [23]. Copulas are 
constructs from probability theory for assessing and sam- 
pling from multivariate distributions. They are widely 
used in finance and ecology [40,41]. For any absolutely 
continuous multivariate cumulative distribution function 
(cdf ) F(x\, . . . , %d) with marginal cdf s F/(#/), i = 1, . . . , d, 
joint density function f(x\, . . . ,Xd) and marginal density 
functions^- (xf), i = 1, . . . , d, we decompose 

...,Xd) = C (Fi(*i), . . . ,Fd(x d ))M x i)'- • -'fd( x d)> 

(14) 

where c (u\, . . . , Ud) is a unique copula density function 
defined on [0,1]^ with uniformly distributed marginals 
on [ 0, 1]. This copula function covers the full dependency 
structure of the variables. In other words, every joint dis- 
tribution function can be decomposed into the marginal 



behavior of its individual variables and a function covering 
its dependency structure [42]. The MH proposal function 
then generates problem specific proposals with an accord- 
ing dependence structure drawn from a pair copula distri- 
bution. Fitting the copula distribution was done in preruns 
containing 1,000,000 unthinned samples each. They were 
generated for each investigation and model separately. 
For back-and-forth conversion of the prerun samples and 
proposals [23], we naturally applied the according prior 
distributions of the models at hand. Choosing different 
conversion functions is possible, but affects the sampling 
performance. Before starting the final MCMC sampling 
procedure, the maximum a posteriori parameter esti- 
mates were computed by simulated annealing and used 
as initial MCMC sampling values. This makes a burn-in 
period dispensable. For thinning the Markov chains, i.e. 
for drawing approximately independent samples in the 
MCMC procedure, we applied the autocorrelation-based 
Effective Sample Size (ESS) proposed by Kass [43]. The 
ESS holds the number of samples left when the Markov 
chain is thinned such that two consecutive samples can be 
considered approximately independent. The copula-based 
MH approach is able to deal with the dependence struc- 
ture in the high dimensional sampling space and allows 
for high proposal acceptance rates at simultaneously high 
ESSs. Finally, all Bayes factors were computed based on 
30,000 proposals of the copula-based MH algorithm at 
each x n throughout all applications. 

Results and discussion 

In this section, we present the results of our analysis. We 
address the question which model is superiorly fitting the 
data. First, several general results, such as investigation 
dependency of the Bayes factor and effects of parameter 
correlations are shown, before turning to the results of the 
model selection, and their consequences for the HMGU 
and ICRP models. 

Investigation specificity of transfer rates 

In radiation protection the transfer rates for the biokinet- 
ics of radionuclides in the human body are derived from 
data collected in various independent experiments [5]. We 
measured plasma and urine levels in 16 different investiga- 
tions. This poses the question whether the models should 
be compared based on the complete dataset, or whether 
statistical evaluation should be done for each investiga- 
tion individually. While the former approach results in 
one overall Bayes factor, the latter yields 16 investiga- 
tion specific, not directly comparable Bayes factors. All 
investigations follow the same pulse-like time courses in 
the transfer compartment y\ while the excretion rates in 
the urine compartment 3/7 exhibit an exponential decay 
behavior (Figure 2). However, zirconium tracer concentra- 
tions showed up to a 50-fold difference between maximal 
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Figure 2 The experimental data. Plasma and urine data for investigations 1 -1 6 on log-log-timescale. 



plasma concentrations, i.e. for investigation 10 (1.616%) 
and 6 (0.033%). 

To test the hypothesis whether the diversity in concen- 
tration also effects transfer rates and therefore the esti- 
mated Bayes factors, we pairwise compared the posterior 
sample marginals of the MCMC run (corresponding to the 
samples of r = 1) for parameter x-j of the ICRP model 
between all investigations by means of a Kolmogorov- 
Smirnov test. Here x-j was chosen as it directly affects 
the observed plasma levels [8]. Except for one pair, all 
p-values were < 6 • 10 -8 , meaning that the chance of 
falsely rejecting the hypothesis of comparable marginals 
is negligible. Therefore, as the posterior marginal distri- 
butions are quite different, it can be deduced that the 
basis for the Bayes factor, the joint posterior distribution, 
can differ quite strongly w.r.t. the individuals. This indi- 
cated that each investigation should be treated separately. 
Nevertheless, in order to infer the transfer rates of an 
average subject, as currently used by the ICRP, the con- 
catenated data had to be used. We therefore compared the 
HMGU and ICRP model based on both the concatenated 
data V = {V\, . . . , V\^} and, in order to account for the 
biological diversity, the individual investigation specific 
datasets V{ (i = 1, . . . , 16). This could also be the basis 
for further analysis of influence factors such as weight 
or gender. 

MCMC-based parameter estimation 

Throughout, the analysis was based on 30,000 propos- 
als for each of the 30 temperature levels in all 17 cases 
(one for each investigation and one for V). For the HMGU 
model the average ESS including one standard error, i.e. 
taken over all temperature levels and investigations, is 
5832 ± 405. In case of the ICRP model we obtained 
in average 5808 zb 252 (approximately independent) 



samples from the Markov chains. This naturally implied 
high acceptance rates for both models. The sampling 
procedure thus captured the power posteriors very 
well. 

From the posterior samples, we could derive new cred- 
ible intervals for the parameters at hand as well as a 
new MAP estimate for an average subject which can 
be used if single parameter values are required (see 
Additional file 1 section 4.1). As can be seen in Figure 3, 
the fit of the time courses covered the data, indicat- 
ing that both models are in principle able to fit the 
data. This justifies our ODE approach with additive 
noise. However, from the fits alone it is not obvious 
which model is superior. Note that the credible inter- 
vals in Figure 3 represent only the uncertainty based on 
the parameters, in contrast to measurement uncertain- 
ties accounted for by the cr^s and cr^s, which are not 
shown. Clearly, this uncertainty in the parameters is spe- 
cific to the individual investigations or the complete data, 
respectively. 

Parameter correlations and model identifiability 

The posterior probabilities of both the HMGU and ICRP 
model showed strong correlation between the parame- 
ters X7 and #8 throughout all investigations. The esti- 
mated Kendalls rs based on the preruns were thmgu — 
0.8027 ± 0.01 and t IC rp = 0.3452 zb 0.02. This can be 
explained as follows: At time point t = 0 the stomach 
compartment 3/9 is the only compartment with non-zero 
Zr concentration. It is exclusively connected to the small 
intestines yio in all models. Therefore, all Zr compounds 
have to pass through jio, which further on distributes 
them to the observed transfer compartment y\ via X7 or 
to the upper large intestines ys via #8- Aberrations in one 
of the parameters X7 or xs thus have a direct effect on 
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time [d] time [d] 

Figure 3 Posterior time courses. Sample median (solid line) and 90% credible interval (CI, shaded area) for the numerical solution of the time 
courses based on the r = 1 HMGU (blue) and ICRP (red) MCMC samples for the complete plasma data (A), urinary excretion rate over time of the 
complete data (B), plasma data of exemplary investigation 15 (C), and urinary excretion rate overtime of exemplary investigation 15 (D) on a log-log 
scale. The median and CI represent the uncertainty in the parameters, in contrast to measurement uncertainty (not shown). Colored markers are the 
data points. The median and the 90% credible interval were computed pointwisely at each time point over all MCMC-based solutions. For 
readability we truncated plasma plots at 1 • 10 _5 [%] and urine plots at 1 • 1 0 _6 [ 



the amount of zirconium predicted for y\ . This affects the 
according posterior distributions. The same effect is found 
for the complete data V (compare pairwise scatterplots in 
Additional file 1 section 4.2). Despite the parameter 
dependencies, the posterior distributions of the ICRP and 
HMGU model are identifiable for all 16 investigations, this 
is, the investigation specific maximum a posteriori esti- 
mates are well denned and inferable (cf. Additional file 1 
section 4.3). 

Bayesian model comparison 

Using the concept of thermodynamic integration we com- 
pared the HMGU and the ICRP model based on (i) the 
concatenated data V = {T>i, . . . ,V\^\ and (ii) the indi- 
vidual investigation specific datasets T>{ (i = 1, . . . , 16). 
This results in a total of 17 Bayes factors. We found that 
all Bayes factors favored the HMGU model; in 14 out of 
the 17 cases even decisively (cf. Table 1, second column, of 
this section and section 4 of Additional file 1). 

In order to take a closer look at the contribution 
of the plasma and urine data to the above results, we 
computed additional Bayes factors based on the likeli- 
hoods £f(a* Jt|£>/) and Jt|X>/) individually. Here, 
i = 1, ... ,16, ALL and k e {I,H}, where / repre- 
sents the ICRP and H the HMGU model. The time 



courses already suggested better coverage of plasma data 
by the HMGU model (Figure 3 above and section 4.4 of 
Additional file 1); for urine the situation is not that clear. 
This was confirmed by the Bayes factors (see Additional 
file 1 section 4 for sampling details): all 17 Bayes fac- 
tors based on plasma data favored the HMGU model; 
in ten cases even decisively (Table 1, third column). For 
the urine data, three investigations slightly favored the 
ICRP model (Table 1, fourth column). In summary, all 
decisive Bayes factors are in favor of the HMGU model. 
While the HMGU model was never decisively rejected, 
the ICRP model was decisively rejected in the major- 
ity of cases. Hence, in depth analysis showed that the 
HMGU model is superior to the ICRP model with respect 
to zirconium processing in the human body. This not 
only holds investigation-specifically, but also based on the 
complete data. We additionally considered an extension 
of the HMGU model (see Additional file 1 section 1.2 
and 4) which, however, did not show any significant 
improvements. 

Differences in radioactive 95 Zr retention in bone predicted 
by the HMGU and ICRP models 

In internal exposure monitoring, biokinetic models were 
used to predict the organ retention or daily excretion 
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Table 1 Bayes factors 



Inv. 


B H,I 


B H,I 


B H,I 


1 


7.17 • 10 1 


7.12 • 10 1 


1.05 


2 


1.15 • 10 2 


2.93 • 1 0 2 


3.94- 10 3 


3 


5.95 • 10 4 


5.23 • 10 4 


1.34 


4 


1.07 • 10 3 


2.64 • 1 0 3 


3.47- 10 1 


5 


2.19- 10 2 


4.73 • 1 0 2 


1.34- 10 2 


6 


4.64- 10 3 


3.93 • 10 3 


2.38- 10 3 


7 


2.18- 10 2 


2.30- 10 2 


1.34- 10 3 


8 


3.75 • 10 1 


1 .28 • 1 0 2 


0.22 


9 


4.62- 10 2 


2.32- 10 2 


0.18 


10 


8.62- 10 2 


1.16- 10 2 


0.20 


11 


1.17 • 10 5 


1.81 • 10 1 


2.94- 10 3 


12 


1.78- 10 2 


5.48 


1.14- 10 1 


13 


7.19- 10 2 


1.41 • 10 1 


4.41 


14 


3.58- 10 1 


7.43 


9.77 


15 


6.29- 10 3 


2.17- 10 1 


1.60- 10 2 


16 


6.22 • 10 2 


1.34- 10 1 


1.20- 10 4 


ALL 


1.20- 10 11 


3.43 • 10 4 


4.73 • 10 7 



Bayes factors for the HMGU versus the ICRP model {B^) for investigation 1 , . . . , 1 6 and the complete data model (ALL) as well as the according Bayes factors for the 
plasma (B P HI ) and urine {B U HI ) data. The HMGU model is favored substantially, when B Hj > 3 and decisively, when B HI > 1 00. Also, }/B HS = B' IH . 



of incorporated radionuclides [44]. With an interval of 
120 days the radioactivity of 95 Zr possibly incorporated 
by occupational workers was routinely monitored by 
whole body counters. Depending on the intake route, the 
radiation dose of bone surfaces or colon was taken as reg- 
ulatory limit for a decision if an individual is requested 
for person-specific monitoring [45]. In this monitoring 
procedure, the biokinetic model structure and parame- 
ters are used implicitly in the background. The organ 
retention function is the solution of the model in each 
compartment; the organ doses are directly related to 
the integral of radioactivity of 95 Zr in source organs 
over 50 years. 

In order to compare the retention of 95 Zr as pre- 
dicted by the ICRP and HMGU models, the 90% credible 
intervals for the time courses in the bone compartments 
were calculated based on the posterior samples. It is 
found that there is a significant difference between the 
two models (Figure 4), where for the ICRP model we 
added the concentrations in the two bone compartments. 
The time courses were derived for stable isotopes of Zr 
and thus had to take the radioactive decay of 95 Zr with 
half-life of 64.032 d [46] into account. The decrease of 
retention in bone using the HMGU model consequently 
reduces the radiation dose estimate in bone in com- 
parison to the ICRP bone dose which is currently used 
in monitoring. 



Retrospective dose assessment 

Internal doses due to incorporated radionuclides have to 
be estimated with the help of biokinetic models based 
on indirect measurements, using for example bioassays 
for blood or urinary excretion. Normally, bioassay or 
in vivo data (e.g. radioactivity accumulated in skull or 
knee detected by a partial body counter) are measured 
after an accidental intake of radionuclides. Uncertain- 
ties of estimated doses are significant and have a large 
impact on the remediation and thus action costs. In con- 
trast to conventional uncertainty analysis as performed 
in [7], our Bayesian approach naturally integrates the 
uncertainties of measured data and parameters simulta- 
neously. This trait of the Bayesian approach is useful as 
it provides an estimate for the intake and its credible 
intervals. 

For example, if the urinary excretion after accidental 
exposure is measured, we can compute credible intervals 
for the initial intake of radionuclide 95 Zr by exploiting 
the posterior distribution together with the linearity of 
the HMGU model. In order to be as general as possi- 
ble we used the posterior samples based on the complete 
data. Given a posterior sample a measurement ytj 
in [jJLg/d] for the urinary excretion rate of zirconium 
at time point t corresponds to a unique solution c^it) 
of the HMGU ODE system. Due to the linearity of the 
ODEs, the initial concentration cy/(0) is by definition 
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10"* 10"' 10 u time[d] 10 ' 10* 10° 

Figure 4 Zirconium retention in bones. Median (solid lines) as well as 90% credible intervals (shaded areas) for the retention of 95 Zr in the bone 
compartment(s) as predicted by the HMGU and ICRP models, taking into account radioactive decay. 



zero for all except the stomach compartment 3/9. The lat- 
ter reads 3/9(0) = y f 7 • 100%/cr^ (£) where c 9 ^{t) denotes 
the value of (£) in the stomach compartment at time 
point t. Now, assuming that for arbitrary posterior sam- 
ples the measurement y f 7 is contained in the 90% 
credible interval of the solution c x H(t) with initial con- 
dition 3/9(0) as given above, lower and upper bounds for 
credible regions of the initial amount of zirconium taken 
in at to = 0 emerge. These are based on the posterior 
samples. The estimated extrapolation factors for multipli- 
cation with a urine measurement (in [/xg/d]) after time 
t (in [h]) are contained in Table 2 and yield the initial 
amount of zirconium contained in the stomach at to = 0. 
For example, if an amount of yif = 50/xg/d was mea- 
sured after two days, we find from Table 2 that the 90% 
credible interval for the ingested amount lies between 
0.029g and 0.059g. Since the above calculations are based 
on non-radioactive Zr isotopes, the results have to be 
corrected for dose assessment with respect to radioac- 
tive decay of the radionuclide in question, i.e. in many 
' 95 Zr. 



cases 



Table 2 Urine predictions for the HMGU model 



Time t 


6h 


12h 


18h 


24h 


30h 


Ibffor IC 


1233.91 


1 820.44 


2614.48 


3369.70 


4100.16 


mffor IC 


1 763.73 


2225.90 


3153.70 


4228.19 


5340.23 


ubfforIC 


2512.54 


2832.49 


3978.27 


5650.86 


7516.00 


Time t 


36h 


42h 


48h 


54h 


60h 


Ibffor IC 


4778.27 


5352.64 


5800.77 


6153.80 


6450.74 


mffor IC 


6364.76 


7250.67 


7977.31 


8557.87 


9006.97 


ubfforIC 


9122.11 


10655.01 


11878.81 


12960.61 


13903.07 



Shown are the lower bound factor (Ibf), median factor (mf), and upper bound 
factor (ubf) for multiplication with a urine measurement (in [/zg/d]) after time f 
(in [h]) on a 60h grid yielding the initial intake concentration (IC) at to = 0. 



Conclusions 

We were the first to evaluate two competing biokinetic 
ODE models for zirconium processing in the human body 
after ingestion. These models were the current model 
recommended by the International Commission on Radi- 
ological Protection (ICRP) and a model developed by 
the Helmholtz Zentrum Munchen (HMGU). The anal- 
ysis was based on a Bayesian approach, i.e. individual 
Bayes factors for 16 investigations as well as a Bayes 
factor based on the concatenated dataset. In order to 
obtain reliable Monte Carlo sampling results, we com- 
bined the numerically stable thermodynamic integration 
with an efficient copula-based Metropolis-Hastings algo- 
rithm. In summary, the HMGU model was unequivocally 
superior with 14 of 17 Bayes factors being even decisive 
when compared to the well-established ICRP model. Also, 
when restricting the data on plasma and urine measure- 
ments only, we found that the HMGU model was clearly 
favored. The HMGU model thus best covers human 
data. 

In contrast to the ICRP model, the HMGU model pre- 
dicted a delayed accumulation of zirconium in bones 
which might be experimentally tested in animals in the 
future. Furthermore, we showed that the HMGU model 
can be applied for retrospective dose assessment, where 
the initially ingested amount of zirconium can be recon- 
structed including credible intervals from ex post urine 
measurements. This provides a simple hands-on tool 
that facilitates the decision if measures have to be taken 
in case of accidental exposure. In future applications 
the superior HMGU model together with its posterior 
samples can readily be used as the basis for dose esti- 
mation in internal dosimetry. The Bayesian framework 
for the compartmental model analysis developed in the 
present work is directly applicable to a personalized 
dose assessment and the uncertainty quantification if a 
person-specific monitoring is requested. More generally, 
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the presented methodology is suitable for any ODE-based 
model selection task, such as the modeling of protein 
signaling, gene regulation, or drug processing [47], nowa- 
days frequently put forward in systems biology [48,49] or 
pharmacogenetics [50]. 



Additional file 



Additional file 1: Supplementary information. Supplementary 
information, including the detailed ODE systems for both models, the prior 
information used for the inference and more detailed evaluation of the 
sampling results, among them additional time course plots for the single 
investigations and scatterplots for the evaluation of parameter 
dependencies. Furthermore, we provide an identifiability analysis for all 
models and a model variant of the HMGU model including its evaluation 
via Bayes factors. 
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